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Abstract 

Molecular dynamics (MD) and continuum simulations are carried out to investigate the influ- 
ence of shear rate and surface roughness on slip flow of a Newtonian fluid. For weak wall-fluid 
interaction energy, the nonlinear shear-rate dependence of the intrinsic slip length in the flow over 
an atomically flat surface is computed by MD simulations. We describe laminar flow away from 
a curved boundary by means of the effective slip length defined with respect to the mean height 
of the surface roughness. Both the magnitude of the effective slip length and the slope of its 
rate-dependence are significantly reduced in the presence of periodic surface roughness. We then 
numerically solve the Navier-Stokes equation for the flow over the rough surface using the rate- 
dependent intrinsic slip length as a local boundary condition. Continuum simulations reproduce 
the behavior of the effective slip length obtained from MD simulations at low shear rates. The slight 
discrepancy between MD and continuum results at high shear rates is explained by examination 
of the local velocity profiles and the pressure distribution along the wavy surface. We found that 
in the region where the curved boundary faces the mainstream flow, the local slip is suppressed 
due to the increase in pressure. The results of the comparative analysis can potentially lead to the 
development of an efficient algorithm for modeling rate-dependent slip flows over rough surfaces. 

PACS numbers: 83.50.Rp , 47.60.Dx , 47.61.-lc , 02.70.Dh , 02.70.Ns , 47.15.-x 
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I. INTRODUCTION 



An adequate description of fluid flow over a solid boundary with surface roughness on 
multiple length scales often requires resolution of fine microscopic details of the flow structure 
while retaining peculiarities of a macroscopic picture The most popular and practically 
important example is the problem of transport through heterogeneous porous media, which 
is an inherently multiscale system [2i|. Modeling of the fluid flow at multiple scales is a 
nontrivial problem itself; however, it becomes much more difficult if one also needs to in- 
corporate slip boundary conditions into equations of motion. Typically, a local slip appears 
if the wall-fluid interaction is sufficiently low, but the value of slip length (an extrapolated 



distance o 
structure 



" the velocity profile to zero) generally depends not only on the surface energy and 



n 



y, 
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11| . Therefore, a correct form of 



but on shear rate as well [5|, 
the Navier-Stokes equation can be obtained only if the combined effect of surface roughness 
and shear rate on the slip length is taken properly into account. It is intuitively clear that 
the surface roughness and shear rate have an opposite impact on the slip length. A precise 
evaluation of relative contributions of these two factors is the subject of the present study. 

The boundary conditions for the flow of monatomic fluids over atomically flat surfaces can 
be described in terms of the intrinsic slip length, which is deflned as the ratio of the shear 
viscosity to the friction coefficient at the liquid/solid interface. The friction is determined by 



the strength of the interaction potential and molecular-sea^ 



wall. Recent MD studies 
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e corrugations of the crystalline 
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2l| have reported values of 



the intrinsic slip length up to about a hundred molecular diameters for weak wall-fluid 
interaction energies. It is well established that the formation of commensurate structures of 
liquid a„d .oM pha.e. at the interface lead, to .tic. bo.„da.. condition. Q.Q. The .Hp 
length decreases with increasing bulk pressure in the flow of nonwetting simple fluids over a 
smooth substrate [l^. The MD simulations have also demonstrated that, depending on the 
wall-fluid interaction energy and the ratio of wall and fluid densities, the slip length is either 
relatively small (below a couple of molecular diameters) and shear-rate- independent jl^ or 
larger than several molecular diameters and increases with shear rate 15|, |20||. Although 



several functional forms for the rate-dependent intrinsic slip length have been suggested 



for monatomic fluids 
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22l |. as well as for polymer melts 
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1261, 12 



the rate-dependent boundary conditions were not used for modeling slip flows in complex 
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geometries. 

In contrast to the description of the flow over smooth boundaries (with microscopic sur- 
face roughness) by the intrinsic shp length, it is more appropriate to characterize the flow 
away from macroscopically rough surfaces by the effective slip length, which is usually de- 
fined with respect to the location of the mean roughness height. Recently, we investigated 



the behavior of the effective slip length in 



for simple fluids 



28l | and polymer melts 



aminar flow over periodically corrugated surfaces 



251 ]. Both MD and continuum simulations have 



shown that the effective slip length decreases with increasing amplitude of the surface rough- 
ness. In the continuum analysis, the local (rate-independent) slip length is modified by the 



presence of boundary curvature [29|, [SOj and, therefore, becomes position-dependent along 



the wavy wall. A direct comparison between MD and continuum simulations at low shear 
rates demonstrated that there is an excellent agreement between the effective slip lengths 
for large wavelengths (above approximately 30 molecular diameters) and small wavenum- 

]. In the present paper, the analysis of the flow over a periodically 



bers ka < 0.3 
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corrugated surface is extended to higher shear rates and to the situation where the local slip 
length at the curved boundary is a function of both the radius of curvature and shear rate. 

In the last decade, a number of studies have focused on the development of the hybrid 
continuum-atomistic methods to simulate fluid flows near smooth solid boundaries 



31 



32|, 



rough walls |33| . superhydrophobic surfaces 3J], corner flows 35|], and flows near a moving 



contact line 



36|, 



37j. In multiscale algorithms, the time and length scales accessible in 



molecular simulations are considerably extended by employing the continuum method in 
the bulk region and applying the MD technique only near the interfaces. The two methods 
are coupled via constrained dynamics in an overlap region by matching mass, momentum and 
energy fluxes or velocity fields. The full MD simulations are then performed to validate the 
hybrid scheme. The drawback of the hybrid methods, however, is the cumbersome procedure 
for coupling of the two methods in the overlap region. The approach pursued in the present 
paper is different from that in the hybrid schemes. We first compute by MD simulations the 
intrinsic slip length as a function of shear rate in the flow over an atomically smooth surface. 
The Navier-Stokes equation is then solved numerically (in the whole computational domain) 
for the flow over a rough surface with the rate-dependent intrinsic slip length used as a local 
boundary condition. We found that the full MD simulations recover the continuum results 
at low and intermediate shear rates and small-scale surface roughness. 
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In this paper, the dynamic behavior of the effective shp length is investigated in steady 
shear flow over a periodically corrugated surface using molecular dynamics and continuum 
simulations. Both methods show that the effective slip length is nearly constant at low shear 
rates and it gradually increases at higher shear rates. The slight discrepancy between the 
effective slip lengths obtained from MD and continuum simulations at high shear rates is 
analyzed by examination of the local velocity and density proflles, as well as pressure and 
temperature distributions along the wavy surface. It is also shown by MD simulations that 
the singular behavior of the intrinsic slip length (in a flow over a smooth wall) at high shear 
rates is suppressed by the surface roughness. 

The rest of this paper is organized as follows. The details of molecular dynamics and 
continuum simulations are described in the next section. The shear rate dependence of the 
intrinsic and effective slip lengths and the results of the comparative analysis are presented 
in Sec. mil The results are briefly summarized in the last section. 



II. THE DETAILS OF NUMERICAL METHODS 
A. Molecular dynamics simulations 

The system geometry and the deflnition of the effective slip length are illustrated in Fig.[Tl 
The total number of fluid monomers in the cell is flxed to Nf = 16 104. The fluid monomers 
interact through the truncated pairwise Lennard- Jones (LJ) potential 



VlAt) = As 



(T\12 /o-\6 



r J \r ) - 



(1) 



where e is the energy scale and a is the length scale of the fluid phase. The cutoff radius for 
the LJ potential is set to = 2.5 cr for both fluid-fluid and wall-fluid interactions. The size 
of the wall atoms is the same as fluid monomers. Wall atoms are flxed at the lattice sites 
and do not interact with each other. 

The motion of the fluid monomers was weakly coupled to a thermal reservoir via the 



Langevin thermostat 38|]. To avoid bias in the shear flow direction, the random force and 
friction terms were added to the equations of motion in the y direction, perpendicular to 
the plane of shear [l^ . The equations of motion for fluid monomers in all three directions 
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are as follows: 



.,. dxi 



myi + niTiji = - + fi , (3) 

-V^ (4) 



where F = 1.0 is a friction coefficient that regulates the rate of heat flux between the fluid 
and the heat bath, and is the random force with zero mean and variance 2mrkBTS{t) 
determined from the fluctuation-dissipation relation. Unless otherwise specifled, the temper- 
ature of the Langevin thermostat is set to T = 1.1 e/A;^, where fc^ is the Boltzmann constant. 
The equations of motion were integrated using the Verlet algorithm 39| with a time step 
At = 0.002 r, where r = ^/mo^Je is the characteristic time of the LJ potential. 

The intrinsic slip length was computed at the lower flat wall composed out of two layers 
of the fee lattice with density Pw= 3.1 a^'^ and nearest-neighbor distance d = 0.77 a between 
atoms in the xy plane. The partial slip is permitted along the lower boundary due to 
relatively low wall-fluid interaction energy ey^f = 0.5e. The no-slip boundary condition is 
enforced at the flat upper wall with a lower density 0.67 cr~^ by allowing the formation of 
commensurate structures of liquid and solid phases at the interface. The nearest-neighbor 
distance between wall atoms is 1.28 o" in the (111) plane of the fee lattice. The wall-fluid 
interaction energy e ensures the no-slip condition at the upper wall holds even at high shear 
rates. The cell dimensions are set to Lx = 42.34 a and Ly = 10.0 a parallel to the conflning 
walls. The MD simulations were performed at a constant fluid density (p = 0.81 a^'^) en- 
semble, i.e., the relative distance between the walls was always flxed to =41.41(7 in the 
z direction. Periodic boundary conditions were applied along the x and y directions. 

The effect of surface roughness on slip flow was investigated in the cell with a lower 
sinusoidal wall with wavelength A = L^. and amplitude a = 1.82 a (see Fig.[T]). Special care 
was taken to make the lower corrugated wall with a uniform density pw= 3.1 a^'^ (by in- 
cluding additional rows of atoms along the y direction) so that the lattice spacing along the 
sinusoidal curve remained (i = 0.77cr. 

The fluid monomer velocities were initialized using the Maxwell-Boltzmann probability 
distribution at the temperature T = 1.1 e/fc^. After an equilibration period of about 2x 10^ r. 
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the velocity and density profiles across the channel were averaged for about 2 x 10^ r within 
bins of Lr X Ly X Az, where A2; = 0.01cr. The location and dimensions of the averaging 
regions for the computation of the local velocity and density profiles will be specified in 
Sec. lITTin 



B. Continuum simulations 



The two-dimensional, steady-state and incompressible Navier-Stokes (NS) equation is 
solved numerically using the finite element method. The equation of motion based on these 
assumptions is written as follows: 

pu ■ Vu = - Vp + yuV^u, (5) 

where u = ui + vj is the velocity vector, p is the density of the fluid, and p and /i are the 
pressure field and fluid viscosity, respectively. 

The penalty formulation was implemented for the incompressible solution to avoid de- 



coupling of the velocity and pressure terms 40|. In the penalty formulation, the continuity 



equation, V ■ u = 0, is replaced with a perturbed equation 

V-u = -|, (6) 

where A is the penalty parameter, which enforces the incompressibility condition. For most 
practical situations, where computation is performed with double-precision 64 bit words, a 



penalty parameter (A) between 10"^ and 10^ is sufficient to conserve the accuracy 40|]. After 
substitution of the pressure term in Eq. into Eq. the modified momentum equation 
can be rewritten as follows: 

pu ■ Vu = AV(V ■ u) + pV^u. (7) 
The Galerkin method with rectangular bilinear isoparametric elements is used to integrate 



the Navier-Stokes equation [40|, |41j. 



The boundary conditions are applied at the inlet, outlet, and upper and lower walls. 
The periodic boundary conditions are implemented along the x direction for the inlet and 
outlet. The boundary condition at the upper wall is always no-slip. A partial slip boundary 
condition is applied at the lower wall. In addition to the global coordinate system {x,z), a 



local coordinate system {t, n) is defined along the lower curved boundary with unit vectors 
t and n representing the tangential and normal directions, respectively. The tangential 
component of the fiuid velocity in the local coordinate system can be computed as 

Ut = Lo[{n-W)ut + Ut/Rl (8) 

where R is the local radius of curvature and Lq is the intrinsic slip length at the flat 
liquid/solid interface [s^. The radius of curvature R is positive and negative for concave 
and convex regions, respectively. The second-order forward-differencing scheme was used to 
compute accurately the normal derivative of the tangential velocity {n ■ V) Ut at the lower 
boundary. 

At the beginning of the simulation procedure, the boundary conditions at the upper and 
lower walls are set to no-slip. Then, the equations of motion are solved implicitly and the 
tangential velocities at the lower boundary are updated according to Eq. ([H]). In the next 
step, the updated velocities at the lower boundary are used as a new boundary condition and 
the equations of motion are solved again. The iterative process is repeated until a desired 
accuracy is achieved. The convergence rate of the solution was controlled by applying an 
under-relaxation factor of 0.001 for the boundary nodes. For all results reported in the 
present paper, the continuum simulations were performed with 150 x 150 grid resolution. 
However, the accuracy of the results was also checked by performing a set of simulations 
using a finer grid 180 x 180. The maximum error in the effective slip length due to the grid 
resolution is Lcfj/L^ = 0.003. The same error bars were reported in the previous study of 
laminar flow over a periodically corrugated surface with either no-slip or (rate- independent) 



slip boundary conditions 



|4l|. 



The averaged error of the solution is deflned as 

r I n" _ h 

i=l I I 

where Np is the number of nodes in the computational domain, u" is the magnitude of 
the velocity at node i and time step n, and u""*"^ is the magnitude of the velocity in the 
next time step. The solution is converged when error ^ 10^^ and the tangential velocities 
along the lower boundary satisfy Ut = Liocai^, where the local slip length is Liocai = 



{L^'-R-')-' m 
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III. RESULTS 



A. MD simulations: the intrinsic and effective slip lengths 

We first consider flow in the cell with flat upper and lower walls. The averaged velocity 
profiles are plotted in Fig. [2] for the selected upper wall speeds. The profiles are linear 
throughout the channel except for a slight curvature near the lower wall at U = 5.0cr/r. 
Note that with increasing upper wall speed, the slip velocity at the lower wall increases 
monotonically up to about 2(j/t at U = 5.0(T/r. The flow velocity at the top boundary 
remains equal to the upper wall speed due to the commensurability of liquid and solid 
structures at the interface and sufficiently high wall-fluid interaction energy. The shear rate 
was estimated from the linear slope of the velocity profiles across the entire cell. 



The fluid viscosity computed from the Kirkwood relation 



4^ is plotted in Fig. [3] as a 



function of shear rate in the cell with flat walls. In agreement with the results of previous 



MD studies on shear flow of simple fluids 15|, l20|, l2l|] , the shear viscosity /i = 2.15±0.15£:rcr 



remains rate-independent up to 7r ~ 0.072, which corresponds to the upper wall speed 
U = 5.375 cr/r. It is also known that at high shear rates the fluid temperature profiles become 
nonuniform across the cell and the heating up is larger near the interfaces |20l . l2ll . |2J] . To 
ensure that the fluid viscosity is not affected by the heating up, we performed two additional 
sets of simulations at the upper wall speeds U = 1.0 a/r and U = 3.5 cr/r and increased the 
temperature of the Langevin thermostat with an increment of 0.05 e/ks- In the range 
of the thermostat temperatures reported in Fig. [31 the shear viscosity is independent of 
temperature for both upper wall speeds. In Sec. lIII C| the results of simulations performed 
at higher temperatures will be used to estimate the effect of pressure on the slip length. 

The intrinsic slip length Lq in the flow over a flat lower wall was determined from the 
linear extrapolation of velocity profiles to zero with respect to the reference plane located 
0.5 cr above the fee lattice plane. The variation of the intrinsic slip length as a function of 
shear rate is presented in Fig. HI At low shear rates 'jt < 0.005, the slip length is nearly 
constant. The leftmost data point reported in Fig.lHis 7r ^ 1.7 x 10~^ at the upper wall 
speed U = 0.01cr/r. At higher shear rates, the intrinsic slip length first gradually increases 
and then grows rapidly as it approaches 7r ~ 0.072 when U = 5.375 a/r. The behavior of 
the slip length in the flow over a flat wall can not be well described by the power law function 
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suggested in Ref. (l5|. Instead, a ninth-order polynomial function was used to fit the data 
(see the blue curve in Fig.H]). In the continuum analysis presented in the next section, the 
polynomial function will be used as a local (rate-dependent) boundary condition, Eq. (jS]), 
for the fiow over a periodically corrugated wall. 

When the upper wall speed becomes greater than [/ = 5.375 cr/r, the rate-dependence of 
the intrinsic slip length passes through a turning point at the maximum shear rate 7r ~ 
0.072, and the slip length increases (Lq ^ L^) while the shear rate decreases, leading to a 
negative slope of the rate-dependent curve (not shown). In this regime the slip velocity at 
the lower fiat wall is greater than the fiuid thermal velocity = ksT /m. The behavior of 
the intrinsic slip length at very large slip velocities was not considered further in the present 
study. 

Next, the infiuence of surface roughness on the effective slip length is investigated in 
the flow over the lower corrugated wall with wavenumber ka = 27ra/A = 0.27. The velocity 
profiles, averaged over the period of corrugation X = Lx, are linear in the bulk of the fllm and 
curved near the lower wall within about 2 a from the bottom of the valley (e.g., see 25l. l41|). 
Therefore, the effective slip length and shear rate were extracted from a linear flt to the 
velocity proflles inside the bulk region 10 ^ z/a ^ 30. The effective slip length as a function 
of shear rate is also plotted in Fig. HI Both the magnitude and slope of the rate-dependence of 
the effective slip length are reduced in comparison with the results for the flat wall. Following 
the rate-independent regime at 7r < 0.01, the effective slip length increases monotonically 
up to a maximum value Leff(0.14r^^) ~ 16.6 cr at U = 8.0a/T. 

At higher shear rates 7r > 0.14, the effective slip length starts to decrease with further 
increase in shear rate (see Fig.H]). We found, however, that the fluid temperature (density) 
becomes anomalously high (low) on the right side of the corrugation peak. For example, 
the fluid temperature in the region 0.5 < x/A < 0.9 varies up to Ths/s ~ 1.6 — 1.85 
at U = 12.0cr/r (the rightmost point in Fig.H]). Also, the amplitude of the characteristic 
density oscillations normal to the curved boundary at x/A ~ 0.43 is reduced by about 40% 
at U = 12.0 a/r in comparison with the fluid layering at equilibrium conditions (not shown). 
In the next section, the comparison between the effective slip lengths estimated from MD 
and continuum simulations is studied at shear rates 7r < 0.14. 
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B. Comparison between MD and continuum simulations 



In this section, we describe the behavior of the rate-dependent effective shp length com- 
puted from the solution of the Navier-Stokes equation. All variables in continuum simu- 
lations are normalized by the LJ length cr, time r, and energy e scales. The system is 
designed to mimic the MD setup. All nondimensional parameters in continuum simulations 
are marked by the (~) sign. The cell dimensions are set to L^ = 42.34 and 17^ = 40.41 in 
the X and z directions, respectively. Similar to the MD setup, the wavenumber of the lower 
corrugated wall is fixed to ka = 0.27. This value was chosen based on the previous analysis 



of the effective slip length in 
rough surface. It was shown 



aminar flow of simple 



281] and polymeric 



25| fluids over a 



28| that at low shear rates there is an excellent agreement 



between the effective slip lengths extracted from MD and continuum simulations when a 
local (rate-independent) intrinsic slip length is used as a boundary condition along the wavy 
wall with A > 30cr and ka < 0.3. 

It is important to emphasize that in continuum simulations the intrinsic slip length in 
Eq. dH]) is a function of the total shear rate at the curved boundary. The total shear rate is 
computed as a sum of the normal derivative of the tangential velocity dut/dn and the ratio 
of the slip velocity to the radius of curvature Ut/R. The effective slip length is extracted 
from a linear part of the velocity profiles in the bulk region 10 ^ 5 ^ 30. 

FigureH] shows the effective slip length extracted from the continuum solution with the 
local rate-dependent intrinsic slip length at the lower corrugated boundary. The upper wall 
speed is varied in the range 0.025 ^ t/ ^ 6. As expected, the effective slip length obtained 
from the continuum simulations agrees well with the MD results at low shear rates 7r < 0.01 
where Lq is nearly constant. With increasing shear rate up to 7r < 0.1, the effective slip 
length increases monotonically. There is an excellent agreement between the two methods 
at intermediate shear rates 0.01 ^ 7t < 0.04. At higher shear rates, 0.04 ^ 7t < 0.1, the 
continuum predictions slightly overestimate the MD results. The maximum discrepancy in 
the effective slip length is about 2.5 a at 7r ^ 0.1, which corresponds to the upper wall 
speed U = 6.0. For U > 6.0, the local shear rate at the crest of the lower boundary, 
[{•n ■ V) "Uj + Ut/R], becomes larger than the highest shear rate, 7r ~ 0.072, achieved in the 
MD simulations for the intrinsic slip length (see Fig.H]), and, therefore, the local boundary 
condition, Eq. ([8]), cannot be applied. 



10 



C. A detailed analysis of the flow near the curved boundary 

In order to determine what factors cause the discrepancy between the effective shp lengths 
obtained from MD and continuum simulations at high shear rates, we performed a detailed 
analysis of the fluid temperature, pressure, density and velocity profiles near the lower bound- 
ary. The fluid temperature in the first fluid layer along the lower wall is plotted in Fig. [5] for 
the indicated upper wall speeds. As expected, at low shear rates •yr < 0.019 {U < l.Ocr/r), 
the fluid temperature remains equal to the value T = 1.1 e/kB set in the Langevin thermo- 
stat. With increasing upper wall speed, the temperature increases and eventually becomes 
nonuniformly distributed along the curved boundary (see Fig.E]). Somewhat surprisingly, 
we flnd that at high upper wall speeds U ^ 5.0 a/r the fluid temperature is lowest above 
the crest even though the slip velocity in this region is expected to be higher than in the 
other locations. 

The normal pressure distribution along the lower corrugated wall is shown in Fig. [6] for 
the selected upper wall speeds. Each data point represents the ratio of the averaged normal 
force on the wall atoms from the fluid monomers (within the cutoff radius) to the surface 
area 0.77 cr x 10.0 cr in the ty plane. The small variation of the pressure at equilibrium is due 
to the boundary curvature, i.e., the number of fluid monomers within the cutoff distance 
from the wall atoms is on average slightly larger above the peak than inside the valley. With 
increasing upper wall speed, the normal pressure increases in the valley {x/X < 0.12 and 
x/A > 0.54) while it becomes smaller than the equilibrium pressure on the right slope of the 
peak (see Fig.E]). The normal pressure distribution extracted from the solution of the NS 
equation (not shown) is qualitatively similar to the proflles presented in Fig.O However, in 
the region x/A > 0.54 the MD data overestimate the continuum predictions for the pressure 
distribution at f/ > 3.0 a/r due to the increase in fluid temperature (see Fig. [5]). It is 
clear from Fig. [6] that at higher upper wall speeds there is a considerable deviation from 
the equilibrium pressure, which can affect the local intrinsic slip length along the curved 
boundary (see discussion below). We flnally comment that the normal pressure on the flat 
lower wall {ka = 0) at equilibrium is equal to P ~ 2.36 e/cx^ (which agrees well with the 
pressure on the locally flat surface at x/ A ^ 0.5 for f/ = in Fig. [6]) and it increases up to 
P 2.52 6 /a^ at the highest upper wall speed U = 5.375 a/r. 

A snapshot of the flow near the lower corrugated wall is shown in Fig. [3, where the 
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rectangular boxes illustrate the location of the averaging regions along the surface. The 
dimensions of the boxes are d = 0.77 a, 5.0 cr, and Ly = 10.0 cr along the curved boundary {xz 
plane), normal to the curved boundary {xz plane), and in the y direction, respectively. The 
velocity and density profiles discussed below were averaged inside the boxes within bins of 
thickness 0.01 cr in the direction normal to the surface. 

In Figure[8]the averaged density profiles are plotted in regions I-IV for U = and 6.0 cr/r. 
In both pronounced fluid layering is developed near the lower wall. In the stationary 

case, Fig.[8](a), the density profiles are almost the same in all four regions. Similar to the 
results of previous MD studies on slip flow over smooth surfaces 20|, l21[, when the upper 
wall speed increases, the contact density (magnitude of the first peak) is reduced due to 
a finite slip velocity along the curved boundary [see Fig.[H](b)]. Among regions I-IV, it is 
expected that the local slip velocity is highest above the crest (region II) and lowest at the 
bottom of the valley (region IV), and, therefore, it is not surprising that the contact density 
in the region IV is larger than in the region II [see Fig.[H](b)]. However, this correlation 
does not hold in regions I and III, where it is expected that at finite Re the slip velocity in 
the region I would be larger than in the region III; but, as shown in Fig.[8](b), the contact 
density in the region I is slightly larger than in the region III. Instead, we correlate this 
behavior with the nonuniform pressure distribution along the curved boundary where the 
local normal pressure in the region I is higher than in the region III (see Fig.|6]). Note also 
that the location of the second and third fluid layers in regions II and III is slightly shifted 
away from the surface and the magnitude of the peaks is reduced due to the lower pressure 
in the region 0.2 < a;/A < 0.5 at f/ = 6.0 cr/r in Fig.O 

The local tangential velocity profiles normal to the curved boundary are presented in 
Fig. [9] for the upper wall speed U = 1.0 a/r. Note that the profiles are almost linear away 
from the surface except for a slight curvature in the MD velocity profiles near the wall. The 
slip velocity above the crest (region II) is larger than in the other regions and Ut at the 
bottom of the valley (region IV) is the lowest, as expected. The tangential velocity on the 
left side of the peak (region I) is larger than on the right side (region III) due to the inertial 
effects 41|. The tangential velocity Ut and its normal derivative dut/dn extracted from MD 
simulations in regions I, III, and IV agree rather well with the continuum results. Despite 
a slight overestimation of the continuum velocity with respect to the MD results above the 
crest (region II), the effective slip length is essentially the same in both methods for the 
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upper wall speed U = 1.0 u/r (see the data at jr ^ 0.019 in Fig.Hj). 

FigureflOl shows local profiles of the velocity component normal to the curved boundary 
for the upper wall speed U = l.Ocr/r. Within statistical uncertainty, the normal velocities 
at the boundary are zero in both MD and continuum simulations. In regions I and IV the 
normal velocity is negative, while in regions II and III the velocity is positive away from the 
boundary. The negative and positive signs of the normal velocity are expected in regions I 
and III because of the positive and negative slopes of the boundary with respect to the shear 
flow direction. Although the boundary slope at the crest and at the bottom of the valley is 
parallel to the upper wall, the local symmetry of the flow with respect to x/X = 0.25 and 
0.75 is broken at finite Reynolds numbers {Re ~ 15.6 for U = 1.0 a/r in Fig. fTUi) . and, as a 
result, the normal velocity in regions II and IV is not zero. Furthermore, the MD velocity 
profiles exhibit an unexpected oscillatory pattern away from the surface. These oscillations 
correlate well with the locations of the minima among the first three fluid layers [e.g., see 
Fig.[H](a)]. Visual inspection of the individual molecule trajectories indicates that the fluid 
monomers mostly undergo thermal motion within the moving layers but occasionally jump 
back and forward between the layers. The net flux between the layers, however, is not zero 
and its sign is determined by the direction of the flow near the boundary. 

The tangential velocity profiles averaged in regions I-IV at a higher upper wall speed 
f/ = 6.0 cr/r are shown in Fig.fTTl The flow pattern near the lower boundary is similar 
to that described in Fig.[9]for U = 1.0 cr/r, but the relative difference between the MD and 
continuum velocities is larger in each region which results in a maximum discrepancy between 
the effective slip lengths of about 2.5 a at 7r ~ 0.1 in Fig. HI We tried to estimate the local 
intrinsic slip length from the MD velocity profiles in regions I-IV and to compare Lq with 
the continuum predictions. However, the normal derivative of the tangential velocity cannot 
be computed accurately because of the curvature in the velocity profiles and ambiguity in 
determining the range for the best linear fit. The resulting error bars associated with the 
linear extrapolation of the MD velocity profiles in Fig. [11] are greater than the maximum 
difference between MD and continuum results for Les at high shear rates in Fig.Hl Instead, 
a more reliable data for the local intrinsic slip length (see Fig. [T2|) are obtained by computing 
the friction coefficient {kf = ^/Lq), which is the ratio of the wall shear stress to the slip 
velocity. The agreement between the MD and continuum results is quite good on the right 
side of the peak in Fig. [121 The intrinsic slip length computed from MD simulations is about 
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3 — 4 0" smaller than the continuum Lq on the left side of the peak due to higher fluid pressure 
in that region (see also discussion below). 

The normal velocity profiles averaged in regions I-IV along the lower boundary are dis- 
played in Fig. [12] for the upper wall speed U = 6.0 cx/t. In each region, the amplitude of 
oscillations in the MD velocity profiles is increased with respect to the mean flow predicted 
by the continuum analysis. Similar to the correlation between density and velocity profiles 
observed in Fig. [TO] for U = l.Ocr/r, the locations of the maxima in the normal velocity pro- 
files in Fig. [13] correspond to the minima in the density profiles shown in Fig.[8](b). At the 
higher upper wall speed, U = 6.0a/T, however, the location of the first peak in the velocity 
profile above the crest (region II) and on the right side of the peak (region III) is slightly 
displaced away from the boundary due to a small shift in the position of the second and 
third fluid layers in regions II and III [see Fig.[8](b)]. 

The effect of pressure variation along the curved boundary is not included in the rate- 
dependent intrinsic slip length, Eq. ([8]), used in continuum simulations. It was shown in 
Fig. [6] that at higher upper wall speeds, U > 3.0 cr/r, the fluid pressure along the lower 
corrugated wall deviates significantly from the equilibrium pressure. Next, we estimate the 
influence of pressure on the intrinsic slip length in shear flow over a flat wall by increasing the 
temperature of the Langevin thermostat. The MD simulations described in Sec. lIII A] were 
repeated at the upper wall speeds U = 1.0 a/r and 3.5 a/r to examine Lq at low {jr ^ 0.018) 
and high {'jr ^ 0.055) shear rates. The results are presented in Fig. [14] for the indicated 
temperatures of the Langevin thermostat. As the bulk pressure increases, the intrinsic slip 
length is reduced by about 3a at U = 3.5 cr/r and about a at U = 1.0 a/r. Note, however, 
that under these flow conditions the shear viscosity remains constant (see Fig. [3]), and the 
variation of the friction coefficient {kf = fJ^/Lo) is determined only by the intrinsic slip 
length. 

The maximum increase in pressure and the difference in Lq reported in Fig. [14] (a) correlate 
well with the variation of pressure in Fig. [6] and the intrinsic slip length in Fig. [12] on the 
left side of the peak at U = 6.0 a/r. We conclude, therefore, that the discrepancy between 
the effective slip lengths computed from MD and continuum simulations at high shear rates 
(shown in Fig.[l]) is caused by the increase in fluid pressure on the left side of the peak 
(where the solid boundary faces the mainstream flow) and, as a result, the local intrinsic 
slip length is reduced. Thus, a more accurate continuum description of the slip flow over 
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a curved boundary can be achieved by including the effect of pressure and temperature on 
the rate-dependent intrinsic shp length. 



IV. CONCLUSIONS 



In this paper, molecular dynamics simulations were used to determine the intrinsic slip 
length as a function of shear rate in steady flow over an atomically fiat surface. It was found, 
in agreement with previous MD studies, that the intrinsic slip length is nearly constant (rate- 
independent) at low shear rates and it grows rapidly as the shear rate approaches a critical 
value. In the presence of periodic surface roughness, the magnitude of the effective slip 
length is significantly reduced and its dependence on shear rate becomes less pronounced. 

For the same surface roughness, the Navier-Stokes equation was solved numerically with 
the rate-dependent boundary conditions specified by the intrinsic slip length obtained from 
MD simulations. The continuum simulations reproduce accurately the MD results for the 
effective slip length at low and intermediate shear rates. The small difference between the 
effective slip lengths at high shear rates was analyzed by performing a detailed analysis of 
the local velocity fields and pressure distribution along the curved boundary. We found 
that the main cause of the discrepancy between MD and continuum results at high shear 
rates is the reduction of the local intrinsic slip length in the region of higher pressure where 
the boundary slope becomes relatively large with respect to the mainstream flow. These 
findings suggest that the continuum analysis of the flow over a rough surface at high shear 
rates should take into account the effect of pressure on the rate-dependent intrinsic slip 
length. 
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FIG. 1: (Color online) The projection of x and z coordinates of the fluid monomers (open circles) 
confined between solid walls (filled circles). The atoms of the lower stationary wall are uniformly 
distributed along the sinusoidal curve with wavelength A, amplitude a, and density p^^, = 3.1 o"^^. 
The shear flow is induced by the flat upper wall with the density 0.67 o""^ moving with a constant 
velocity U in the x direction. The effective slip length Leg is computed by extrapolation of the 
velocity profile to Ux = 0. 
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FIG. 2: (Color online) Averaged velocity profiles for the tabulated upper wall speeds in the cell 
with flat upper and lower walls. The solid lines are the best linear fits to the data. The vertical 
axes indicate the location of the fee lattice planes. 
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FIG. 3: (Color online) Shear viscosity in the cell with flat upper and lower walls when the tem- 
perature of the Langevin thermostat is set to T = l.le/kB (□)• Fluid viscosity as a function of 
the thermostat temperature 1.1 ^ Tks/s ^ 1-4 for the upper wall speed U = 1.0(t/t (A) and 
1.1 ^ TkB/e ^ 1.35 for U = 3.5 a/r (0). 
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FIG. 4: (Color online) The intrinsic slip length L^/a as a function of shear rate for atomically 
flat walls (□). The blue solid curve is a ninth-order polynomial fit to the data. The effective slip 
length Leflf/o" as a function of shear rate for the flow over the corrugated wall with wavenumber 
ka = 0.27 (o). The effective slip length computed from the solution of the Navier-Stokes equation 
with the local rate-dependent slip length (•). 
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FIG. 5: (Color online) Temperature of the first fluid layer T{x) (in units of e/ks) along the lower 
corrugated wall {ka = 0.27) for the indicated upper wall speeds. The vertical dashed lines denote 
the location of the crest (x/X = 0.25) and the bottom of the valley (x/A = 0.75). 
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FIG. 6: (Color online) The normal pressure Pn (in units oi e/a"^) along the lower corrugated wall 
{ka = 0.27) for the selected upper wall speeds. The vertical dashed lines indicate the location of 
the crest {x/X = 0.25) and the bottom of the valley {x/X = 0.75). 
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FIG. 7: (Color online) A snapshot of the fluid monomers (open circles) near the lower corrugated 
wall (filled circles) with wavenumber ka = 0.27. The rectangular boxes denote averaging regions 
on the left side of the peak (I) (x/A ~ 0.05), above the crest (II) (x/A ~ 0.25), on the right side of 
the peak (III) (x/X ~ 0.45), and at the bottom of the valley (IV) (x/A ~ 0.75). The unit vectors 
t and n indicate the tangential and normal directions at the boundary. 
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FIG. 8: (Color online) Fluid density profiles averaged inside regions HV (shown in Fig. [7]) near 
the corrugated lower wall with wavenumber ka = 0.27. The upper wall speeds are U = (a) and 
U = 6.0a/T (b). 
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FIG. 9: (Color online) Averaged tangential velocity profiles inside regions HV obtained from MD 
simulations (open symbols) and the solution of the Navier-Stokes equation (filled symbols). The 
location of the averaging regions near the lower corrugated wall is shown in Fig. [71 The upper wall 
speed is U = 1.0ct/t. 
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FIG. 10: (Color online) Averaged normal velocity profiles inside regions I-IV (shown in Fig.[7I) for 
the upper wall speed U = l.Oa/r. The profiles are extracted from MD simulations (open symbols) 
and the solution of the Navier-Stokes equation (filled symbols). The horizontal dashed lines indicate 
the location of the minima in density profiles (see text for details). 
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FIG. 11: (Color online) Tangential velocity profiles averaged inside regions I-IV (shown in Fig.[7I) 
near the lower corrugated wall. The upper wall speed is U = 6.0(t/t. The velocity profiles are 
extracted from MD simulations (open symbols) and the solution of the Navier-Stokes equation 
(filled symbols). 
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FIG. 12: (Color online) The intrinsic slip length along the lower corrugated wall {ka = 0.27) for 
the upper wall speed U = 6.0 a/r. The data are obtained from MD simulations (open symbols) 
and the solution of the Navier-Stokes equation (filled symbols). The vertical dashed lines denote 
the position of the crest (x/A = 0.25) and the bottom of the valley (x/A = 0.75). 
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FIG. 13: (Color online) Averaged normal velocity profiles inside regions I-IV obtained from MD 
simulations (open symbols) and the solution of the Navier-Stokes equation (filled symbols). The 
upper wall speed is C/ = 6.0(T/r. The horizontal arrows indicate the location of the minima in 
density profiles between the first and second fluid layers in regions I and III. 
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FIG. 14: (Color online) The intrinsic slip length Lo/a as a function of the bulk pressure P (in 
units of s/a^) in the cell with flat upper and lower walls. The temperature of the Langevin 
thermostat is varied in the range 1.1 ^ T/c^/e ^ 1.35 for the upper wall speed U = 3.5a/T (a) and 
1.1 TkB/e ^ 1.4 for C/= l.Ocr/r (b). 
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